#!/bin/csh -f
#goal: read *FreqHca.bin   
#      then draw HCA-height picture
#  by: wangh 2015/3/18
#-----------------------------------------
#set debug
set cmd = $0
set cmd = $cmd:t
#-----------------------------------------
if ( $#argv == 0 ) goto usage

switch ( $#argv )
   case 2:
   breaksw
  default:
usage:
  echo " "
  echo "Usage : $cmd EndTime(yyyymmdd) ExpNam(FR or WR )"
  echo " "
  exit(2)
  breaksw
endsw
  set enddate_UTC = $1
  set ExpNam = $2
#----------------------------------------
#Time information
#----------------------------------------
  echo "---------------------------------------------"
  echo "UTC:  ${enddate_UTC}"
  echo "---------------------------------------------"
#-----------------------------------------
#set homedir = "/public/home/wangh29/data2019-2021/code/StatisPol/"
set homedir = "./data"
set datadir = "${homedir}/${ExpNam}/${enddate_UTC}FreqHca.bin"
set HgtFile = "${homedir}/${ExpNam}/${enddate_UTC}PolQuan.bin"
set outdir  = "./svg/fig9.${ExpNam}.${enddate_UTC}HCA-hgeight${ExpNam}"
#-----------------------------------------
#Parameter set 
#-----------------------------------------
set type = "x11" #pdf,x11,ps,ncgm
set FontS = 0.025
#-----------------------------------------
#Create plot.ncl
#-----------------------------------------
  cat >obs_radar.ncl << EOF
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;  These files still have to be loaded manually
;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
  begin
;-------------------------------------------
   wks_type ="${type}"
   if("${type}".eq."png")then
    wks_type@wkWidth = 2500
    wks_type@wkHeight = 2500
   else
    wks_type@wkPaperHeightF  =11 
    wks_type@wkPaperWidthF = 11 
   end if 
   wks  = gsn_open_wks(wks_type,"${outdir}")
   cmap = RGBtoCmap("~/script/colorbar/radar_HCA.rgb")
   gsn_define_colormap(wks,cmap)
;-------------------------------------------------
  colors = ispan(1,10,1)
  labels  = (/"GC","BS","DS","WS","CR","GR","BD","RA","HR","RH"/)
;-------------------------------------------
; read hgt data
  setfileoption ("bin", "ReadByteOrder", "BigEndian")
  xyz  = fbinrecread("${HgtFile}",0,(/3/),"integer") 
  nz = xyz(2)
  Hgt1d = fbinrecread("${HgtFile}",45,(/nz/),"float") 
  hgt1d = Hgt1d/1000. ; m to km
;------------------------------
; read HCA frequency 
  xy  = fbinrecread("${datadir}",0,(/2/),"integer") 
  nz  = xy(1)
  nh  = xy(0)
  print(xy)
  hca1d = fbinrecread("${datadir}",1,(/nh/),"float") 
  temp  = fbinrecread("${datadir}",2,(/nz,nh/),"integer")
  AvgFrlv = fbinrecread("${datadir}",3,(/1/),"float")
  Avg20  = fbinrecread("${datadir}",4,(/1/),"float")
  AvgFrlv = AvgFrlv/1000.
  Avg20 = Avg20/1000.

  freqT = new((/nz,nh/),"float",-999.)
  freqH = new((/nz,nh/),"float",-999.)
  freqT = 0.
 ; print(temp(13,:))
  do k=0,nz-1
   total = sum(temp(k,:))
   if(total.gt.1)then
    freqT(k,:) = int2flt(temp(k,:))/int2flt(total)*100.   
   end if 
   do n=0,nh-1
     print(n+" hgt "+sprintf("%6.2f",hgt1d(k))+" "+labels(n)+" "+sprintf("%6.2f",freqT(k,n)))
     freqH(k,n) = sum(freqT(k,0:n))
   end do 
  end do 
  freqH!0   = "hgt"
  freqH&hgt = hgt1d
  ;print(freqH(13,:))

;----------------------------------
; Set some basic resources
  res   = True
  res@gsnDraw    = False
  res@gsnFrame   = False
;  res@gsnMaximize              = True    ; Maximize plot in frame
;  res@vpWidthF                 = 0.3     ; Make long and
;  res@vpHeightF                = 0.9     ; narrow
;---Set axes limits. Add extra space for X max.
  res@trXMinF                  = 0.0
  res@trXMaxF                  = 100.
  res@trYMinF                  = 1.
  res@trYMaxF                  = 15.
;---Put city labels on Y axis
  res@tmYRLabelFontHeightF  = ${FontS}
  res@tmYLLabelFontHeightF  = ${FontS}
  res@tmXBLabelFontHeightF  = ${FontS}
  res@tiYAxisString            = "Height (km)"
  res@tiXAxisString            = "Frequency (%)"
  res@gsnLeftString  = " a) ${ExpNam}"
  res@gsnLeftStringOrthogonalPosF =-0.12 
  res@gsnLeftStringFontHeightF  = ${FontS}
  res@tmXMajorGrid                = True    ; Turn on grid lines
  res@tmXMajorGridLineDashPattern = 2       ; Dashed lines
  res@tmXMajorGridThicknessF      = 1.0
  res@tmYMajorGrid                = True    ; Turn on grid lines
  res@tmYMajorGridLineDashPattern = 2       ; Dashed lines
  res@tmYMajorGridThicknessF      = 1.0


;---Create base plot with just outlined bars
  base_plot = gsn_csm_blank_plot(wks,res) 
  delete([/res@tmXMajorGridLineDashPattern,res@tmXMajorGridThicknessF/])
  plotB = new((/nz,nh/),"graphic")
;---Create plot of shorter bars
  bres           = True
  bres@gsEdgesOn = True      ; Outline the polygons (bars)
 do k=1,nz-2
  HgtS = (hgt1d(k)+hgt1d(k-1))/2.
  HgtE = (hgt1d(k)+hgt1d(k+1))/2.
  ybar = (/HgtS,HgtS,HgtE,HgtE,HgtS/)
  do n=0,nh-1
    if(n.eq.0)then
     xbar  = (/0.,freqH(k,n),freqH(k,n),0.,0./)
    else
     xbar  = (/freqH(k,n-1),freqH(k,n),freqH(k,n),freqH(k,n-1),freqH(k,n-1)/)
    end if 
    bres@gsFillColor         = colors(n)
    bres@gsEdgeThicknessF    = 0.0
    bres@gsEdgeColor         = bres@gsFillColor
    plotB(k,n) = gsn_add_polygon(wks,base_plot,xbar,ybar,bres)
  end do
 end do
;--------------------------------------------------------
;------------------------------------
;line resourse
  resp                  = True                ; polyline mods desired
  resp@gsLineColor      = "black"             ; color of lines
  resp@gsLineThicknessF = 2.0                 ; thickness of lines
  resp@gsLineDashPattern= 0
  xpts = (/-100., 100./)
  ypts = (/AvgFrlv,AvgFrlv/)
  dum  =  gsn_add_polyline(wks,base_plot,xpts,ypts,resp)
  ypts = (/Avg20,Avg20/)
  resp@gsLineDashPattern= 1
  dum2  =  gsn_add_polyline(wks,base_plot,xpts,ypts,resp)
  txres = True
  txres@txFontHeightF = ${FontS} 
  txres@txFontColor  = resp@gsLineColor
  text0= gsn_add_text(wks,base_plot,"0~S~o~N~C",10.,AvgFrlv+0.5,txres)
  text2= gsn_add_text(wks,base_plot,"-20~S~o~N~C",10.,Avg20+0.5,txres)

;**********************************************************
; add labelbar to plot
  getvalues base_plot         ; get plot size for use in creating labelbar
  "vpXF"      : vpx
  "vpYF"      : vpy
  "vpHeightF" : vph
  "vpWidthF"  : vpw
  end getvalues

  lbw    = 0.4 * vph           ; Make labelbar size a fraction of the plot.
  lbh    = 0.5 * vpw
  nboxes = nh;dimsizes(res@gsnXYBarChartColors)

  lbres                    = True          ; labelbar only resources
  lbres@vpWidthF           = 0.18 * vph     ; labelbar width
  lbres@vpHeightF          = 0.8 * vpw     ; labelbar height
  lbres@lbBoxMajorExtentF  = 0.8          ; puts space between color boxes
  lbres@lbFillColors       = colors; res@gsnXYBarChartColors ; labelbar colors
  lbres@lbMonoFillPattern  = True          ; Solid fill pattern
  lbres@lbLabelFontHeightF = 0.025         ; font height. default is small
  lbres@lbLabelJust        = "CenterLeft"  ; left justify labels

  gsn_labelbar_ndc(wks,nboxes,labels,0.82,0.7,lbres)

  draw(base_plot)             ; this draws all three XY plots.
  frame(wks)

;-----------------------------------------------

end
EOF

ncl obs_radar.ncl 
rm  obs_radar.ncl

